Entanglement Dynamics in Harmonic Oscillator Chains 
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We study the long-time evolution of the bipartite entanglement in translationally invariant ID 
harmonic lattice systems. We show that for a wide class of Hamiltonians and generic initial states 
there exists a lower bound for the von Neumann entropy which increases linearly in time. This 
implies that the dynamics of harmonic lattice systems can in general not efficiently be simulated by 
algorithms based on matrix-product decompositions of the quantum state. 
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Recently, so-called matrix product states (MPS) have 
received much interest for the numerical simulation of 
one-dimensional quantum many body systems [1], 2\. 
This is because the ground state of fermionic and bosonic 
lattice systems with finite-range interactions and an ex- 
citation gap usually implies an area law of entanglement 
Q , stating that the von Neumann entropy of a partition 
scales with the surface size. In ID the surface area is in- 
dependent on the size of the system resulting in a weakly 
entangled ground state, which can thus faithfully be rep- 
resented by MPS. Even when the excitation gap vanishes, 
i.e. for critical systems, there is only a correction which 
is at most logarithmic in the system size. The situation 
is however quite different for non-equilibrium problems 
as here not only the scaling with size but also with time 
is relevant. With respect to the latter only an upper 
bound derived by Lieb and Robinson exists |4j], which 
states that the von Neumann entropy increases at most 
linear in time. Being an upper bound, it does of course 
not allow to draw any conclusion about the approxima- 
bility of the long-time dynamics of quantum many-body 
systems by MPS. However, it is very often found, that 
the bi-partite entropy does indeed scale linear in time, 
which implies that the required computational resources 
increase exponentially in time. For example it has been 
shown for the spin-i XY model that the entropy grows 
linearly with time after a global quench |5| . On the other 
hand, it was found recently for the case of free fermions 
that the entropy can grow only logarithmically in time 
@, showing that for certain initial states the long-time 
dynamics is accessible with MPS based methods. If for 
free fermions the scaling of the entanglement entropy in 
time is only moderate for certain initial conditions, what 
kind of scaling with time can we expect for the entangle- 
ment entropy of free bosons? 

In the present paper we study the time evolution of the 
entropy in ID bosonic systems that evolve under trans- 
lationally invariant quadratic Hamiltonians with local or 
finite-range couplings. In order to separate the prob- 
lem of size-scaling from the scaling with time, which is 
the subject of interest in this paper, we consider a spe- 
cific class of ID systems where the bi-partite entropy be- 



comes independent on system size. To this end we choose 
as initial state of the time evolution the ground state <£> 
of some local, gapped Hamiltonian Hq. As explained in 
the following, one can intuitively expect that under such 
initial conditions the bi-partite entanglement entropy of 
the time-evolved state will be independent on system 
size. Since for the ground state of local Hamiltonians 
the presence of an excitation gap is sufficient for an area 
law of entanglement [3|, the entropy of the initial state 
is size independent. It is easy to see that the state time- 
evolved under a Hamiltonian H, vE'(i) = exp{— iHt}$, 
is the ground state of a time-dependent Hamiltonian 
H'[t] = cxp{-iHt}H cx_-p{iHt}. The spectrum of H'[t] 
is identical to that of Hq, i.e. it too has an excitation 
gap. Furthermore the Lieb-Robinson bounds guarantee 
that for any fixed time t its coupling matrix elements 
between sites i and j are exponentially small beyond a 
certain distance l c , i.e. for \i — j\ > l c . Thus H'[t] is also 
of finite range Q • As a consequence we can expect that 
the entanglement entropy of the time-evolved state will 
saturate with increasing system size. However it is not 
possible to draw any conclusion about the time-scaling of 
entanglement beyond the limits set by the Lieb-Robinson 
upper bounds. 

In the present paper we show that in contrast to free 
fermions the entropy of the time evolved quantum state 
always grows linearly in time, making a long-time simu- 
lation of ID bosonic systems with MPS based methods 
impossible. 

Translationally invariant harmonic oscillators: To be 
specific we consider a one-dimensional system of ./V 
bosonic oscillators described by N pairs of canonical op- 
erators x = (xi, X2, ...,xn) and p = (pi,p2, ■ ■■,Pn)- The 
oscillators are coupled by a quadratic Hamiltonian of the 
form 
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(x| V |x> 



(1) 



where V is a real, symmetric, positive definite, time- 
independent matrix. We assume translational invariance, 
implying that V is a Toeplitz matrix. Furthermore we 
consider periodic boundary conditions, such that V is 
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circulant. Circulant matrices form a commutative alge- 
bra. Moreover the elements of a circulant matrix can be 
generated from the spectral function X(9), i.e. 



V k i = 



i r 27r 



i(k-l)9 



(2) 



We now want to determine the scaling of the entan- 
glement entropy with time. To this end we have to find 
the time evolution under the local Hamiltonian (JT]). The 
system is assumed to start its evolution at t = from a 
Gaussian state, i.e. 



<I>(x) = m i1 pxi>( --(x\B\x) 



a 
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where B is a real, symmetric, and positive definite 
Toeplitz matrix. Periodic boundary conditions imply 
that B is also a circulant matrix with spectral function 
/3(6). As the initial state is the ground state of a gapped, 
local Hamiltonian, (3(9) is non-zero and regular corre- 
sponding to a non-critical state. Since the Hamiltonian 
of the system is quadratic, the time-evolved state remains 
Gaussian and we have to search for a solution in the form 



*(x,t) = 



1 



N dot A~ 



i/4 exp(--(x| J 4(t)|x) 



(4) 

Here and in the following a tilde denotes the real part, 
i.e. X = x+ 2 ■ By taking into account the symmetry of 
B and after simple calculations one can easily find that 
A (t) obeys the Riccati equation 
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V, 



A{0) = B. 



Its solution can be written as 

1/2 cos (tV 1 / 2 ) B + iV 1 ' 2 sin (tV 1 ^ 2 ) 



A{t) = V L 



(ty 1 /2)yi/2 + isin (tvV 2 )B 



(5) 
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which is again a circulant matrix. The spectral function 
A(9, t) of its real part A can easily be obtained from X(9) 
and P(0): 



A(0,t) 



X (9) cos 2 {tX 1 ' 2 (9)) + f3 2 (9) sin 2 (tX 1 / 2 (9)) ' 

(7) 

Note that if B = V 1 / 2 , the spectral function and thus the 
matrix A(t) becomes time- independent, as in this case 
the initial state is the ground state of the full Hamilto- 
nian. 

Reduced density matrix: Having the solution of the 
Schrodinger equation we can now calculate the reduced 
density matrix of a block of N — n oscillators. The cal- 
culations can be done by partitioning the symmetric ma- 
trices A (t) and A^ 1 (t) into blocks 



A(t) = 



T C 
C T R 



A' 1 (*) 



Q D 
D T P 



(8) 



where T is an n x n and R an (TV — n) ~x(N — n) matrix. 
Similar calculations have been done in |8j and [9| for the 
ground state of a chain of oscillators. After a lengthy 
but straightforward calculation we find for the matrix 
elements of the reduced density operator 



p R (x, x') = Afexp 
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, (9) 



where x = (a; n +i, ....xn), x' = (x' n+1 , ....x' N ) are the co- 
ordinates of the remaining N — n oscillators, 

R C T f- x C C T T- X C* 
r =2 — A = ^I— ' 



and M = ( det P 



1/2 



/(•7r) 2 is a normalization with 
R-C T f- 1 C, (10) 



Here and in the following Y 1 denotes (y^J 

Purity and lower bound for the Entropy: We proceed 
by analyzing the dynamical behavior of the bi-partite 
entanglement. There are several measures of entangle- 
ment between parties of a closed system, examples be- 
ing the von Neumann entropy S = — tr (p^ mp^) and 
the purity trp^ , where the following inequality holds 
S > — lntr[pjj]. It should be noted that — lntrp|j rep- 
resents also a lower bound to all Renyi entropies S a = 
In tr [p R ] with a < 1 as S a > Si = S. 
In order to derive a lower bound for the entropy we 
calculate the purity of ((9]). 

tr [Pr] = J dxdx Vfl (x, x') p R (x', x) . (11) 

The Gaussian nature of (j9]) allows to calculate this inte- 
gral in a straight-forward way: 



tr [Pr] 
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After simple algebra one obtains 
tr [ P r] = 



det [P [R + Z T T- 1 Z 



A 



-1/2 



1/2- 



(12) 
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where Z = (C — C*)/2i. The last inequality follows 
from the fact that Z T T^ 1 Z is a positive definite ma- 
trix. With this we find the following lower bound to the 
von-Neumann entropy 



S > - In det ( / ' A' 



(13) 



In order to facilitate analytical calculations of determi- 
nants, we consider the limits N ^> 1 and N > n ^> l.It 
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can then be shown [10( that in this limit the elements 
of matrices R and P can be generated from the spectral 
functions A(6,t) and A -1 (9,t) respectively. As A(6,t) 
is a regular function i.e. A (9, t) > for any t, we may 
apply the strong Szego theorem [Tljj to calculate the de- 
terminants. According to this theorem 



Cfc 



(14) 



fc=i 



where the c k are Fourier coefficients of In A 1 (8,t), i.e. 
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Cfe = — / d6> In A" 1 
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?,f) exp (— z0fc) . (15) FIG. 1: (Color online) Numerical plot for the sum y^k 
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If ft (9) and A((9) are constant function i.e. the oscillators 
are uncoupled, all Fourier coefhcents (I15|) vanish except 
for Co. In this case (jT4"]) reduces to the trivial bound S > 
( entanglement is never generated). In what follows we 
will consider only the non trivial case when X(9) is a not 
constant. 

In Fig. [1] we have plotted the right hand side of eq.(fT4|) 
numerically evaluated for an initial state with spectral 
function ft(9) = 1 and a Hamiltonian H with spectral 
function X(9) = (c — cos(6*)) 2 . If c > 1, H has a finite 
excitation gap as there are no real zeroth of X(9). If the 
gap vanishes, i.e. for c < 1 the ground state of H becomes 
critical. One clearly recognizes a linear increase with time 
in all cases. That the presence of an excitation gap is 
irrelevant here is not surprising, as the initial state has 
a finite overlap with excited states except in the trivial 
case where it coincides with the ground state. 

For short times one should expect that the sum grows 
quadratically in time The spectral function A (9,t) 
([7]) for small t scales approximately quadratic in t 



A(0,t)^ft(9) [l-(ft 2 (9)-X(9))t 2 



(16) 



The correction to the initial spectral function is propor- 
tional to difference ft 2 (9) — X(9), as was expected. The 
Fourier coefficients c k (p"5|) can then easily be calculated 
Cfe ~ £fc + t 2 Sk ,where and 6k are some constant num- 
bers. From this one can calculate the sum (fT4|) for short 
time which yields 



S > H\ + H^t 2 



(17) 



In the following we will derive an analytic estimate for 
the lower bound to the entropy for large times. Note 
that ft 2 (9) — X(9) corresponds to an initial state that is 
an eigenstate of H and thus has no time evolution at all. 
In any real system, the number of oscillators in the chain 
is finite and therefore to neglect boundary effects in the 
thermodynamic limit it is necessary consider time inter- 
vals t < L/v, where v is the speed for excitations after a 
quench, the so-called Lieb-Robinson speed Q (see also: 



function of time for a Hamiltonian H with spectral function 
A (9) = (c — cos9) 2 and for initial Gaussian state ft (9) = 
1 . The top-most curve (blue) correpsonds to the critical 
Hamiltonian with c = 0.5, the middle curve (magenta) to the 
critical Hamiltonian with c = 1, and the lowest curve (yellow) 
to a gapped Hamiltonian with c = 1.5. One clearly recognizes 
a linear increase with time. The insert shows the quadratic 
short-time evolution. 



[7|), and L is the system size. For the sake of simplicity of 
the derivations we consider a Hamiltonian H with a finite 
excitation gap. The derivation for a non-gapped Hamil- 
tonian is more involved and will not be presented here. 
As noted above the presence of a gap is however irrele- 
vant. In the following we use an alternative expression 
00 2 

for ^2 k \ c k\ which is very useful for numerical and ana- 
fe=i 

lytical calculations. By Parseval's theorem this sum can 



be rewritten as ^~~^fc |c fc | z = 
fe=l 
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Making use of the inequality 
< |x| , \y\ < M one finds 



In 
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M 



s > 



fe=l 



(18) 



where M = max A (8, t), and 



b k = — d9A- 1 (9,t)exp(-i9k). (19) 
27r Jo 

The coefficients have a simple physical meaning: They 
determine the correlations in momentum space over a 
distance k, i.e. {t)\pip i+ k\^ {t)) ~ bfc- With this we 
have 



oo 



(20) 



fe=i 



where we have decomposed A 1 (9,t) in a time indepen- 
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dent and a time dependent term 







(21) 



1 f ,„A(0) ~ P 2 (0) 





;(2U 1/2 (<9)) cos(fc<9). 

(22) 



The term proportional to c; 2 in eq. ([2H|) does not de- 
pend on time and can be disregarded. The sec- 
ond term can be rewritten using the triangular in- 



equality and Parseval theorem to give 



J2 k^kHk (*) 
k=i 



< 



fc=i k=i 
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The time de- 



fe=i 



pendence of this term is thus given by the square root of 
the term containing /ij,(i) 2 . 

We now show that the term ~ /ifc(^) 2 in ea.(|20p is 
bounded from below by a function linear in t. To this 
end we evaluate the integral in eq. ([22l for Hk(t) by the 
method of stationary phase (the role of the large param- 
eter is played by t). The stationary points of the phase 
2A 1 / 2 (0) ± \Q are the solution of 



1 d\(6,t) k 
X 1 l 2 {9,t) d9 ± 1 



0. 



(23) 



As the interaction matrix V is of finite range, the spec- 
tral function A (9) is a trigonometric polynomial of fi- 
nite degree K. By the theorem of Bernstein [ll[ one has 

< K max A {&) , and therefore if k > k max = 
23l) has no real solutions. For 



d\(9) 
d8 



K max X(0) , , 

; v ' t = v a t eq 

^/min MO) 9 ^ 

these values of k all fik (t) are exponentially small in 
agreement with the finite Lieb-Robinson speed. On the 
other hand, when t and k are large but k < v g t one finds 



1 W 

v m— 1 



l) cos ( tGm (l 



(24) 



where F m (-|j and G m (|) are some "nice' 

'k\ „„^^„„+;™„l f AO \\ (o\ a2 



functions. 

F m (|) is proportional to J d6» [A (6») - /3 2 (6>)] which 
quantifies the difference between initial state and ground 
state of H. If both agree, i.e. if j3(9) 2 = A(6»), the 
coefficient vanishes. The integer W is the number of 
stationary points, which must be finite because A (9) 
is a trigonometric polynomial of finite degree. ip m = 
±2 depending of the sign of the second derivative of 

oo 

A 1 / 2 (9) at the stationary points. Thus ^ fc/i 2 (t) > 



k=l 



tVg 

EknUt) 
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By 



highly oscillating terms we arrive at 
w 



t f 



k=l 

Thus 
with 



m— 1 



5 > Ct, 



C 



2M 2 



dxx 



E( F «^))' 



(25) 



(26) 



being a finite time- independent constant. Eq. (|25p is the 
main result of our paper. It constitutes a lower bound to 
the scaling of the entanglement entropy S a , a < 1 with 
time in a one-dimensional system of coupled harmonic 
oscillators. Eq. (|25|) implies that the bond dimension of 
the matrices used in an MPS representation needs to in- 
crease exponentially in time to allow for a faithful rep- 
resentation of the dynamical many body wavefunction. 
This means that in contrast to fermionic systems, where 
at least for certain initial conditions a simulation of the 
long-time dynamics is possible, for harmonic oscillator 
systems this is in general impossible. The discussion can 
be extended to d dimensions. In higher dimensions the 
form of the reduced density matrix is the same as Eq. (j9|) . 
To calculate determinants of Toeplitz matrices one can 
apply the cZ-dimensional Szego theorem (l3| . 
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